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Abstract 

In this work, we present a symplectic integration scheme to numerically 
compute space debris motion. Such an integrator is particularly suitable to 
obtain reliable trajectories of objects lying on high orbits, especially geosta- 
tionary ones. Indeed, it has already been demonstrated that such objects 
could stay there for hundreds of years. Our model takes into account the 
Earth's gravitational potential, luni-solar and planetary gravitational per- 
turbations and direct solar radiation pressure. Based on the analysis of the 
energy conservation and on a comparison with a high order non- symplectic 
integrator, we show that our algorithm allows us to use large time steps 
and keep accurate results. We also propose an innovative method to model 
Earth's shadow crossings by means of a smooth shadow function. In the 
particular framework of symplectic integration, such a function needs to be 
included analytically in the equations of motion in order to prevent numerical 
drifts of the energy. For the sake of completeness, both cylindrical shadows 
and penumbra transitions models are considered. We show that both mod- 
els are not equivalent and that big discrepancies actually appear between 
associated orbits, especially for high area-to-mass ratios. 

Keywords: Space debris; Symplectic integration ; Solar radiation pressure; 
Shadowing effects ; High area-to-mass ratios 
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1. Introduction 

Several works have already shown the great importance of the direct so- 
lar radiation pressure effects on space debris motion, especially in the case 
of high area-to-mass ratios (Liou & Weaver, 2005; Chao, 2005; Anselmo & 
Pardini, 2005; Chao, 2006; Valk et al, 2008; Lemaitre et al., 2009; Valk et 
al., 2009). Moreover, some authors have pointed out that Earth's shadow- 
ing effects should not be neglected and that these phenomena could lead to 
significant short-periodic perturbations (see Kozai 1961, Ferraz-Mello 1972, 
Aksnes 1976 and Valk & Lemaitre 2008). Recent observations (Schildknecht 
et al., 2010) have made it possible to identify the area-to-mass ratio of 274 
uncorrelated objects. The resulting AIUB^/ESA catalogue contains a sig- 
nificant population of objects with area-to-mass ratios larger than 1 m^/kg 
and as high as 86.7 m^/kg. In this framework, perturbations due to solar 
radiation pressure, coupled to Earth's shadows, become essential. This paper 
will focus on the numerical propagation of space debris motion, considering 
both cylindrical shadow and penumbra transition models. 

A second aspect of the space debris dynamics that we decided to take into 
account in this analysis is their possible long lifetimes. If it is obvious, for low 
orbits, that the efficiency of the drag forces cleans the area in few years, it 
is slightly different for higher orbits, especially geostationary ones, on which 
the objects can stay for hundreds of years (see Chao 2005, Valk et al. 2008, 
etc). Up to now the numerical simulations concerning the debris have been 
performed using a wide variety of numerical integrators for short timescales 
and the choice of a symplectic (or quasi-symplectic) integration scheme to 
compute the orbit of space debris has rarely been done (see Breiter et al. 2005 
for an application of a fourth order symplectic integrator of the Wisdom- 
Holman type). Nevertheless, it is well known that these algorithms enable 
the use of large time steps, show excellent energy preservation properties on 
long time scales and are less time-consuming than non-symplectic schemes. 

The introduction of Earth's shadow in this context is a challenging prob- 
lem. Indeed, the actual direct solar radiation pressure reaching the satellite 
has to be computed with a sufficiently smooth function and included directly 
in the equations of motion in order not to introduce numerical errors in the 
symplectic scheme. 

The idea of the introduction of a continuous shadow function equal to one 
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in direct sunlight and zero otherwise has been first proposed in Ferraz-Mello 
(1964) and Ferraz-Mello (1965). This function depends on the angle formed 
by the geocentric Cartesian position of space debris and the dark pole of 
the Earth's terminator. Another approximation of this shadow function has 
been proposed later on in Lala & Sehnal (1969). Then, a way of comput- 
ing shadow boundaries has been proposed in Escobal (1976). It has to be 
noted that, with such approximations, only cylindrical shaped Earth's shad- 
ows could be modeled. Hence, these solutions have been improved further in 
order to take into account penumbra transitions. For example, umbra and 
penumbra cone boundaries have been computed in Escobal (1976). Other de- 
tailed studies exist. For example, penumbra transitions and several physical 
processes in the atmosphere have been added in Vokrouhlicky et al. (1993) 
to provide a realistic shadow crossing model. This theory being really time 
consuming, an approximate version has been presented in Vokrouhlicky et 
al. (1994a). Penumbra phenomena induced by the solar radiation pressure 
from the Earth-reflected sunlight have also been studied in Vokrouhlicky 
et al. (1994b), leading to the conclusion that these effects were much less 
important than the direct solar radiation pressure. A further numerical in- 
vestigation (Vokrouhlicky et al., 1996) has also shown that the oblateness 
of the Earth did not bring significant differences, compared to a spherical 
Earth. In light of this, the oblateness of the Earth has not been included in 
our shadow models. Eventually, we can find another geometrical model of 
the penumbra transition in Montenbruck & Gill (2005), where the compu- 
tation of the degree of occultation of the Sun by the Earth is performed by 
means of apparent radii in geocentric Cartesian coordinates. 

In this paper, we present an innovative theory to model both cylindrical 
and conical shadows (umbra-penumbra transitions) crossings. Our model 
being totally included in the equations of motion of the space debris, this is 
well adapted to our symplectic integration scheme. Besides, our formulation 
has the advantage to build a smooth shadow function. Hence, it can also 
be included in the variational equations whose solutions are required by a 
variety of chaos indicators. It should be noted that it cannot be achieved 
with the formulation from Montenbruck & Gill (2005). 

The paper is organized as follows. In Sec. 2, we describe the forces con- 
sidered in our Hamiltonian model. Sec. 3 is devoted to the presentation 
of the symplectic integrator (Sec. 3.1), its implementation (Sec. 3.2) and 
shows a comparison with an efficient non-symplectic integrator (Sec. 3.3). 
Our shadow modelling theory is then developed in Sec. 4. First (Sec. 4.1), 
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we present the basis of the model for a cyhndrical-shaped Earth's shadow. 
Second, details are given about the adaptation of this theory in the case of 
penumbra transitions in Sec. 4.2. For the sake of clarity, full mathematical 
developments have been moved to different appendices. Eventually, we com- 
pare our penumbra transition model to the one proposed in Montenbruck & 
Gill (2005) in the case of a typical shadow crossing occuring on a geostation- 
ary orbit (Sec. 4.3). We also show the difference between the evolution of 
the orbital elements of a space debris crossing either cylindrical or conical 
Earth's shadows. 



2. Model 

Let us consider the following autonomous Hamiltonian function 
•H(v, A, r, 9) = H^^^^xiy, r) + H.oM + H^^o^otij, 6) 

+'H3body(l") + "H srp V""- / 

where r := [x, y, z) and v are respectively the Cartesian geocentric coordi- 
nates and velocities of the satellite, 9 is the Greenwich sidereal time and A is 
its associated momentum. Our model takes into account the attraction of the 
Earth as a point mass central body ("Hkcpi)? the rotation of the Earth around 
itself ("Hrot) and the perturbations due to the Earth gravity field ("Hgeopot), 
third bodies (mainly the Sun and the Moon) ("Hsbody) and the solar radiation 
pressure ("Hsrp)- 

The attraction of the Earth as a central body is accounted for as 

nj _ 

/tkepl — 

2 r 

where r := ||r||, v := ||v|| and /i = ^M® is the standard gravitational 
coefficient. 

The rotation of the Earth around its axis of smallest inertia is modeled as 

Hrot = OA. 

The complete Earth's potential expressed in the frame rotating around 
the Earth's axis of smallest inertia and with the same angular speed can be 
written as 

f n=Om=0 \ ^ J 

X {Cnm COS m\ + Snm sin m\) 
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where (r, A, 0) are the geocentric spherical coordinates of the sateUite, 
represents the equatorial radius of the Earth, Vnm are Legendre functions 
and Cnm and Snm are the spherical harmonics coefficients. 

Thanks to formulas presented in Cunningham (1970) and later in Mon- 
tenbruck & Gill (2005), we are able to write t/geopot as a function of r 



® n=0 m=0 



Let us remark that Vnm and Wnm depend directly on Cartesian coordinates 
and are defined recursively. Also note that t/geopot is still written in the same 
rotating frame as the one described above. The following simple change of 
variable 

X i-^ X cos 6 + y sin 6 
y I— 7- —X sin 6 + y cos 6 

lets us express f/geopot in the fixed inertial geocentric frame. 

Also note that the zero degree is already present in "Hkcpi- Moreover, 
given that the center of mass coincides with the origin of the reference frame 
Cio = Cu = Sio = Sii = 0. Hence, the sum in t/geopot should start at n = 2. 
It follows that, 



^gcopot(l^) — ^ ^ ^ ^ CnmVnmi^j ^) + Snm^nmi^y 

W n=2 m=0 



The EGM96 model (Lemoine et al., 1970) of degree and order 360 is used as 
Earth gravity model. 

Perturbations due to third bodies (the Sun, the Moon and planets of the 
solar system) are introduced in the Hamiltonian function as 

"Hsbody = — V"" /^j ( 71 n ~ 71 — fT^ ) 

^ \\\T-Vi\\ r,/ •^y 
J \ II II II II / 

where /Xj = QMi and is the geocentric Cartesian coordinates of any third 
body of mass Mj. Position and velocity of the latter are computed by means 
of Jet Propulsion Laboratory (JPL) DE^OS planetary and lunar ephemeris 
(Standish, 1998). This external contribution introduces a time dependence 
in the Hamiltonian, formulated as a quasi periodic function. For simpler 
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models of the solar orbit, limited to three frequencies for example, we could 
have introduced three new variables (linear functions of time) and associated 
conjugated momenta (as it has been done for the rotation of the Earth). 
However it would have considerably increased the number of the differential 
equations to deal with, and the model would stay approximate. This is 
the reason for which we have chosen to work on with a quasi-symplectic 
formulation in this context. 

Eventually, our model also includes the direct solar radiation pressure. 
Three standard physical models are used: the absorption, the reflection and 
the diffusion. Hence, the force induced by the radiation pressure is obtained 
by adding together the elementary forces accounting for each effect. The 
model assumes that space debris are spherical objects and that there is no 
radiation from the surface of the Earth. Detailed information about the 
model can be found in Milani & Gronchi (2009). It leads to the following 
potential: 



where Tq is the geocentric Cartesian position of the Sun, Pr = 4.56 x 
10~^ iV/m^ is the radiation pressure for an object located at a distance of 
1 AU from the Sun, A/m is the area-to-mass ratio of the space debris and 
Qq is equal to the mean distance between the Sun and the Earth (i.e. a© = 1 
AU). The construction of this potential is also explained in e.g. Montenbruck 
& Gill (2005). Let us remark, that, making some assumptions, "Hsrp is writ- 
ten under its conservative form. Moreover, the shadow of the Sun by the 
Earth is not yet considered. 

3. Numerical integration scheme 

3.1. Symplectic integrator 

Let 'H(p, q) = A + eB he an autonomous Hamiltonian with degrees 
of freedom where p, q G are respectively momenta and variable vectors. 
The Hamiltonian vector field can be written as 



•srp 



r — rQll m 




N 



dU ax 



&H ax 

dqj dpj 



X = 



where 
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The solution of this differential equation is given by 

x(t)=e*^«x(to). (1) 

The Campbell-Baker-Hausdorff theorem (Bourbaki, 1972) ensures that we 
can find a general integrator with n steps of the form 

where r is the time span covered by the integration process. Coefficients 
Ci and di are to be chosen carefully, so that the order of the integrator is 
improved on. In this case, integrating H at order m means that we exactly 
evaluate e'^^'^ where 

}C = A + eB + Oir""). 

In Laskar & Robutel (2001), four classes of symmetric symplectic integra- 
tors are presented. 

SABA2n(r) = e^i^^^e'^i^^-s... 
SABA2„+i(r) = e^i^^^e'^i^^-B... 

SBAB2n(r) = e'^i^-^-sgCa^-f^^e'^^rL.B... 

gCn+iTLyigdn+irLgggCn+irL^ 



SBAB2„+i(r) = e'^i^^-8e^2rL^__, 



It is shown that, if e is small enough, one can find a suitable set of coeffi- 
cients {ci,di\ such that 

/C = ^ + ei3 + 0(r2"e + rV). (2) 

As a matter of fact, the order of the error does not only depend on the time 
step anymore but also on e. It yields excellent energy preservation properties 
and enables us to use bigger time steps than with any other non-symplectic 
integration scheme. 
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3.2. Implementation 

As explained in Section 3.1, the Hamiltonian function % has to be spht 
into two integrable parts in order to be expressed as the perturbation of an 
integrable one. Moreover, the S-part must always be smaller than the ^-part 
(i.e. e := is small enough). A convenient possibility is the following 

one: 

?{(v,A,r,^) = ^(v,r,A) + i3(r,^) 

where 



^(v,r,A) 



'Hkcpi(v, r) 



Kot(A) 

+ 'H3body(l") 



^srp(r) 



In this case, each perturbation is put in the i3-part and is of lower order than 
the energy associated to the Keplerian problem (i.e. the ^-part). 

Each time that the operator e'^^'^^^ is applied on the state vector (v, 6', r, A)-^(to) 
the resulting vector corresponds to the solution of the Keplerian problem at 
time epoch to + c^r. This computation is performed analytically in order to 
reduce numerical cost and keep a high accuracy. Moreover, the sidereal time 
is increased linearly 

e{to + qt) = e{to) + Cije 

where 9 is assumed to be constant. 

On the other hand, the solution of e'^''^^=^(v, ^, r, A)^(to) is computed nu- 
merically as the solution of the set of differential equations 



-rfieVr-B(r, 
dieVAB{r,i 









It follows that 

/ v\ 

A 



(to + diT) 



A 



(to 



/ v(to)-rfi£rVri3(r(to),^(to)) \ 
k{ta) - dierV eB{v{t^),e{t^)) 

r(to) 
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3.3. Numerical simulation 

This section aims to prove the numerical efficiency of our integrator. First, 
we show in Fig. 1 the maximum relative error in energy and CPU time 
required by our symplectic scheme to propagate the orbit of a space debris on 
a long time span (500 years), depending on several time steps and integrator 
orders. Initial conditions are a = 42164.140 km, e = 0.1, i = 0.1 rad, 
Q = u = M = rad and the initial Julian date is 2455194.5 days. The model 
includes the geopotential up to degree and order 4. Time steps for these 
simulations have been set to 24, 12, 4, 1 and 1/2 hours. Several observations 
can be made from Fig. 1. On the one hand, the system being fully symplectic 
in this case, the energy is accurately preserved even with large time steps. Let 
us point out that we could not go below the accuracy treshold of 10~^^. On 
the other hand, CPU times are really small, keeping in mind that the orbit 
has been propagated over 500 years. Computations have been performed on 
a E5440 Intel Xeon CPU (2.83 GHz) with 6144 KB cache size. 

The introduction of the gravitational perturbation of the Sun does not 
prevent the energy to be preserved on long time scales. However, the quasi- 
symplecticity of our integrator means that relative variations of the energy 
are about 10~^ for the third body and solar radiation pressure perturbations, 
which is the amplitude of the quasi periodic solar motion. In light of this, 
we compare the relative energy error associated to two different orbits. The 
model still includes the geopotential up to degree and order 4 and we add the 
solar graviational perturbation. The difference appears in the computation 
of the Cartesian position of the Sun. On the one hand, we use the JPL 
ephemeris and, on the other hand, we consider a Keplerian solar motion 
with eccentricity and inclination equal to zero (part of the ephemeris module 
of the NIMASTEP software described below). This lets us reduce the quasi 
periodic motion of the Sun. Results are illustrated in Fig. 2, where the 
SABA4^ integrator has been used. Obviously, a circular and coplanar solar 
orbit reduces the relative error in energy to order 10~^. Fig. 2 also emphasizes 
the excellent energy preservation over the years, even when JPL ephemeris 
are used. 

It has to be noted that the introduction of the JPL ephemeris greatly 
increases the computation cost. As a point of comparison, the CPU time 
needed for this simulation after 500 years with a time step of 4 hours and 
the SABA4 integrator is 22 seconds. It corresponds roughly to twice the 
amount of CPU time corresponding to the same propagation without the 
computation of the Sun's position. 
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Timestep [hour] 



Figure 1: Maximum relative errors in energy (top panel) and CPU times (bottom panel) 
for different integrator orders, as a function of the time step. Initial conditions are a = 
42164.140 km, e = 0.1, i = 0.1 rad, il = uj = M = rad and the initial Juhan date 
is 2455194.5 days. The model includes the geopotential up to degree and order 4. The 
integration has been performed on a time span of 500 years. 
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Figure 2: Comparison of the relative errors in energy as a function of time. Both real and 
simplified (i.e. e© = and i© = rad) orbits of the Sun are considered. Initial conditions 
are a = 42164.140 km, e = 0.1, i = 0.1 radians, 17 = cj = M = radians and the initial 
Julian date is 698382.5 days. The model includes the geopotential up to degree and order 
4 and the graviational perturbation of the Sun. The simulation has been run with the 
SABA4 integrator. 
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Step size [day] 


Step size [s 


CPU time [s 


1/200 


432 


142.01 


1/100 


864 


79.31 


1/86 


1004.65 


70.73 


1/75 


1152 


66.83 



Table 1: CPU times required by NIMASTEP (with ABMIO) with respect to time steps. 

Eventually, we compare our integration scheme to the NIMASTEP soft- 
ware (Delsate & Compere, 2011). NIMASTEP (Numerical Integration of 
the Motion of Artificial Satellites orbiting a TEUuric Planet) is an exten- 
sive tool that allows to integrate numerically the osculating motion of an 
arbitrary object (natural or artificial satellite, space debris...) orbiting a 
central body ((dwarf-) planets or asteroids of the Solar System) taking into 
account a large number of forces, integrators and options. This software has 
been successfully validated and compared to other external softwares. In 
the following comparisons, NIMASTEP is used with the Adams-Bashforth- 
Moulton integrator of order ten (Hairer et al., 1993), hereafter referred to 
as ABMIO. Starting from the same initial conditions, both integrators have 
been used to numerically propagate the orbit. For the sake of completeness, 
all kind of perturbations have been considered. Hence, we take into account 
the geopotential up to degree and order 4, luni-solar perturbations and the 
solar radiation pressure with A/m = 0.01 m^/kg. In this case, the fourth 
order SABA integrator has been used with time steps equal to 4 hours. The 
CPU time required to propagate the initial conditions to 190 years with this 
method is 11.54 s. Table 1 lists each CPU time required by NIMASTEP 
for different time steps. In order to get an idea of the number of steps per 
revolution (about one day for these initial conditions), time steps are given 
in terms of fractions of one day and in seconds. This table let us point out 
that, even if NIMASTEP performs the numerical integration quickly, the 
fact that our symplectic scheme can use large (ten times bigger) time steps 
make it obviously a faster algorithm. 

Then, given that both orbits obtained with our symplectic scheme and 
NIMASTEP are very close to each other, we have decided to show the 
absolute difference between each Keplerian element and for each time step. 
Results can be seen in Fig. 3. For each time step and orbital element, 
the absolute difference has been computed each day. Both softwares having 
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been developed using different units of time and distance, residual errors 
could have been introduced artificially. Anyway, we are pretty confident 
that our comparison does not suffer too much from this kind of numerical 
errors. From Fig. 3, it can be seen that absolute errors are quite small for 
the eccentricity, inclination, longitude of ascending node and argument of 
pericenter, no matter which time step is considered. The situation is slighlty 
different for the semi-major axis. In this case, the error is clearly bigger for 
larger time steps and increases linearly with time. This can be explained 
by the fact that NIMASTEP is a non-symplectic integration scheme and 
that the energy is not fully preserved. Hence, while both integrators turn 
out to be really efficient on short time scales, our symplectic integrator does 
not suffer from a slight drift on the semi-major axis on long time scales. 
Nevertheless, let us remark that NIMASTEP with time steps equal to 432 
and 864 s yields excellent results, the drift on the semi-major axis being 
reasonably low after 190 years (8.4443 x 10"^ km and 2.9404 x 10"^ km 
respectively with times steps 432 and 864 s). One of the main advantages 
of our integrator is that it lets us obtain accurate results even with big time 
steps. As a comparison point, it turns out that the NIMASTEP software 
cannot be used with time steps larger than 1160 s for this particular set of 
initial conditions and perturbations. 

4. Earth's shadow modelling 

The Earth's shadow can first be modeled as a simple cylinder. In this 
case, the Sun is assumed to be far enough from the Earth and the solar 
rays are then supposed to be parallel. The geometry of this problem is 
illustrated in Fig. 4 (top). However, a more realistic model includes the 
penumbra transition which lets us model partial eclipses (see Fig. 4, bottom). 
The distance to the Sun and diameters of both Earth and Sun have to be 
considered to compute the amount of sunlight actually reaching the space 
debris. In the following, we describe an existing method to model cylindrical 
shadows and present our solution (Section 4.1). Then, the same work is done 
in the case of penumbra transition models (Section 4.2). Finally, in Sec. 4.3, 
we compare our model to the solution of Montenbruck & Gill (2005) and 
we present an analysis of differences between orbits obtained with both our 
cylindrical and conical shadow models. 
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Figure 3: Comparison of the absolute difference between each Keplerian elements of the 
orbit obtained by our scheme and by the NIMASTEP software (with ABMIO). Several 
time steps have been used with NIMASTEP: 432 s (green), 864 s (red), 1004.65 s (blue) 
and 1152 s (black). Initial conditions are a = 42164.140 km, e = 0.1 , i = 0.1 rad, 
il = uj = M = rad. The model includes the gcopotcntial up to degree and order 4, solar 
radiation pression (with A/m = 0.01 m^/kg) and luni-solar perturbations. The initial 
Julian date is 2455194.5 days. , . 




Figure 4: Top panel: cylindrieal Earth shadow with solar rays assumed to be parallel when 
reaching the Earth. Bottom panel: Umbra-penumbra model ineluding partial eclipses. 
Angles a and /3 give the geometric difference between the cylindrical and the conical 
models. 
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4.1. Cylindrical shadow models 

In Escobal (1976), a method is proposed to find the orbital entrance and 
exit of a satellite from the shadow of the Earth. More precisely, two non- 
spurious roots of a quartic polynom in the cosine of the true anomaly corre- 
spond to the shadow entrance and exit (more details are given in Appendix 
A). Then the effect of the solar radiation pressure can be switched on and 
off, depending on the angular position of debris on its orbit. Unfortunately, 
it turns out that this shadow crossing model can not be put directly in the 
equations of motion, yielding numerical errors on short time scales. 

In order to avoid numerical errors in the integration process, it is necessary 
to find a smooth function z/(r) equal to one when debris are in direct sunlight 
and zero otherwise. Then, Vr'Hsrp(r) can be multiplied by i^(r) in equations 
of motion so that each shadow crossing is taken into account. 

Starting from a preliminary relation in the model of Escobal (1976), it 
turns out that the satellite is situated in the cylindrical shadow of the Earth 
when 



where the constant 7 has to be fixed according to the required precision. The 
shape of this function is shown for different values of 7 in Fig. 5 (top). It can 
be seen that, the bigger 7, the sharper the function Uc- As a matter of fact, 
a perfect cylindrical shadow model would require 7 to be infinite. However, 
considering double precision floating point standard. Fig. 5 (bottom) shows 
that taking 7 = 10^ is sufficient to represent cylindrical shadow crossings. 
Indeed, the absolute difference between 1 and the function Uc with Sc(r) = 
10~^ is already of order 10~^. 

As explained above, replacing Vr'Hsrp(r) by 1^0(1") Vr'Hsrp(r) in equations of 
motion turns out to be an efficient way to use the fully symplectic integration 
scheme and consider cylindrical-shaped shadows of the Earth. In this case, 
a special attention has to be paid to the integration time step. The latter 
must be small enough to perform some steps inside the umbra zone, which 
only represents a small part of the total revolution time. As illustrated in 
Fig. 7, the cylindrical shadow on a geostationary orbit only lasts around half 





We introduce a new shadow function defined as 



i^c(r) = - |l + tanh[7Sc(r)]| = 



1 in cylindrical umbra 
otherwise 
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Figure 5: Top panel : evolution of the function Vc for different values of Sc and of the 
parameter 7. Bottom panel : absolute difference between one and the function Vc evaluated 
at Sc(r) = 10^^ for different values of the parameter 7. 
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an hour. While this drawback cannot be neglected, it also lets us to use low 
order symplectic integrators. In particular, the SBAB2 integrator is used 
when the shadowing effects are enabled, still keeping highly accurate results. 

4.2. Conical shadow models 

Several attempts have been made to model the penumbra transition. A 
solution is proposed in Escobal (1976) to add umbra-penumbra corrections. 
Basically, one ends up with a more accurate but also trickier shadow function 
which still induces numerical errors in our symplectic integration scheme 
because the radiation pressure does not include Earth's shadow as a smooth 
function. 

Another kind of shadow crossing model can be found in Montenbruck & 
Gill (2005). In this coefficient pm corresponds to the fraction of 

sunlight reaching the debris, based on the angular separation and diameters 
of the Sun and the Earth. Hence, vm is equal respectively to zero and 
one when the debris is in direct sunlight and umbra and corresponds to 
the remaining fraction of sunlight in the penumbra transition phase. It has 
to be noted that this fraction can only be defined in the penumbra cone. 
Hence, it is not possible to handle the function um with a single formula 
at each time and it cannot be used directly within our symplectic scheme. 
Moreover, any stability study requiring the computation of the deviation 
vectors could not be used with this method, um being a piecewise-defined 
function. Nevertheless, vm is kept back as a comparison criterion for our 
further developments. 

In the following, we present an innovative way of modelling umbra and 
penumbra cones crossings during the numerical integration of space debris 
orbit. First, we introduce angles a and /3 representing the difference between 
the umbra cylinder and respectively the umbra and penumbra cones (see Fig. 
4) 



a 



at an ■ 



Rq — Rq 



and (3 = atan ■ 



Rg) + Rg 



with Rq the radius of the Sun. Extending relation (3), it follows that space 
debris are in the umbra cone when 



s„ r 



+ cos a 



— R%, cos^ a + i?m sin a 



< 
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Figure 6: Evolution of the functions Sc, Su and Sp during a shadow crossing on a geo- 
stationary orbit. The penumbra cone is crossed when Sp is negative and the space debris 
lays in the umbra cone while Su is negative. The time spent in the penumbra transition is 
noted At and the difference between Su and Sp at the entrance of the cylindrical shadow 
is denoted by A/i. 



and in the penumbra cone when 



^0 



cos/3 



cos^ (3 — i?0 sin (3 



< 0. 



An example of the evolution of functions Sc, and Sp depending on time is 
shown in Fig. 6. 

Now, we will show that the function can be adapted to include the 
penumbra transition. Actually, the parameter 7 will not be constant anymore 
but will be chosen so that the new shadow function Up is equal to one in direct 
sunlight, starts to decrease in the penumbra cone and is equal to zero in the 
umbra cone. The value of 1 — z/p when the penumbra cone is crossed has to 
be fixed to attain a given precision treshold, denoted a. Hence, we define the 
constant 

S := atanh (1 — a). 

In this paper, 6 is set equal to 8, meaning that the precision treshold a ~ 
2.25 X 10"^. 
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Then, assuming that the time spent in the penumbra transition, At, is 
known, 7 is set equal to 5/ At. With this configuration. 

If. - r 5 

2 



1 + tanh 



At 



(4) 



is smooth and such that 



= 1 
< 1 



1 
a 
a 
a 




if 
if 
if 
if 
if 



^^p(r; 



> 
= 

< 
= 

< 



and 



> 



The main difficuhy hes in the way to estimate At. As a matter of fact, 
this quantity cannot be computed exphcitely before each shadow crossing. 
However, we will show that it can be replaced by a quantity depending only 
on the position of the space debris. Both entrance and exit times spent in 
the penumbra cone being computed exactly in the same fashion, we will only 
explain our method in the entrance case. 

In following developments, each function Sc, s„ and Sp will be expressed as 
functions of the angle between r and Tq. As a first step, we assume that r 
does not depend on 0, meaning that the orbit of the space debris is circular. 
It yields 



Sc{(f>) = rcosi 



■Su(0) = r cos + cos a 



"5p(0) = rcos0 + cos/3 



Ra 



— _R|> cos^ a 



Rm sin a 



R^ cos^ /3 — i?© sin /3 



(5) 
(6) 



(7) 



Let us also define 0i, 02 and 03 respectively as the value of at the entrance 
the cylindrical shadow, umbra and penumbra cones. Hence, the following 
relations hold: 

Sc(0l) = S,(02) = Sp(03) = 0. (8) 

The difference (A0) between 02 and 03 will help us to characterize At. Tech- 
nical details about the computation of A0 are given in Appendix B. It turns 
out that it can be approximated by 



A0 := 



2p — . 

r 



(9) 
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where p := r /tq. Moreover, the hnk between A0 and At can be expressed 

as 



27r 
day" 



-At. 



Hence, 



A0 

At ~ day ~ day. 

2tx Txr 



(10) 



The final step is to show that the difference between and Sp at the entrance 
of the cyhndrical shadow, denoted by A/i, can be used instead of At. Simple 
calculations let us write down A/i as 



R 

Ah = 2pR(s — 



Full details about this relation are to be found in Appendix C. 
In conclusion, we have shown that At can be approximated as 



(12) 



At 



Ah 



From a geometrical point of view, our approximation means that the slope 
of each curve Su, Sc and Sp is close to —1 at the shadow entrance and to 1 at 
the shadow exit with our particular choice of units. 

The same mathematical development can be achieved in the case of non- 
circular space debris orbits (r depends on the angle (p). The calculations are 
presented in Appendix D. 

Going back to equation (4), it is easily seen that we now have 



6271 Rq 
_Ah{r) 



-Sr\r 



(13) 



Let us remark that, in practice. Ah is computed with actual values of a 
and f3 : 

Ah{r) -- 



cosa(r) \/r'^ — -R| cos^ a(r) + R@ sina(r) 



— cos/3(r) a/t"^ — -Rl cos^ /3(r) — R(^ sin /3(r) 
Also note that, with a = /3 = 0, z/p is equal to v^. 
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Figure 7: Typical shadow crossing on a geostationary orbit. This figure shows the com- 
parison between three shadow crossing models, as a function of time. 

4-3. Numerical comparisons 

First, we compare our shadow functions to the one proposed in Mon- 
tenbruck & Gill (2005). It is interesting to study the shape of these functions 
on a single shadow crossing. Considering a geostationary orbit, we show the 
evolution of u^, Up and uu during a typical shadow crossing in Fig. 7. Even 
if the shapes of Up and um are different, both share some common properties. 
It can be seen that both shadow functions cross the cylindrical shadow limit 
at = Up = uyi — 0.5, start to decrease when entering the penumbra cone, 
are equal to zero in the umbra cone and increase again in the penumbra 
exit transition. A numerical comparison based on several shadow crossings 
is proposed later on. 

Short-periodic effects of Earth's cylindrical shadows on the orbital ele- 
ments of space debris have already been studied in Valk & Lemaitre (2008). 
For example, the evolution of the semi-major axis and eccentricity for space 
debris disturbed by the solar radiation pressure and an area-to-mass ratio 
equal to 5 m^/kg is described in Valk & Lemaitre (2008) (Fig. 3). The same 
simulation has been performed in Fig 8 with our conical shadow model. At 
first glance, it is not possible to detect significant discrepancies. 

Nevertheless, both shadow models clearly lead to significantly different 
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Figure 8: Evolution of Kcplcrian elements of the orbit of space debris subject to the 
Earth's central attraction and solar radiation pressure. The initial semi-major axis is set 
at 42164 km, the other elements are set equal to zero and the area-to-mass ratio is equal 
to 5 m^/kg. The shadow function fp is used to model the Earth's shadowing effects. The 
results are in agreement with what is proposed in Fig. 3 in Valk & Lemaitre (2008). 
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debris trajectories. As illustrated in Fig. 9, the absolute difference between 
orbital elements obtained with both models is zero before the first shadow 
season and starts to increase after this period of time. A shadow season 
appears each time that the Sun moves through the orbital plane of motion, 
leading to a succession of shadow crossings. More information about such 
phenomena can be found in Valk & Lemaitre (2008). The area-to-mass ratio 
is equal to 20 m^/kg for this example. After 15 years, the difference between 
both trajectories keeps on increasing and reaches high values, especially for 
the semi-major axis (see Fig. 10, blue curve). In light of this, it turns out that 
cylindrical shadow models are not reliable approximations of conical Earth's 
shadows, especially in the case of space debris with high area-to-mass ratios. 
The altitude of space debris trajectory must also be considered, the higher 
the orbit, the larger the time spent in the penumbra transition. 

In Fig. 10, the results obtained with both shadow functions z/p and z/m are 
represented by the red curve. Given that uu cannot be used with our sym- 
plectic scheme, it has been included in NIMASTEP in order to perform the 
comparison between both methods. Differences between Keplerian elements 
in this case are clearly smaller than in the previous comparison involving Uc 
and Up (Fig. 10, blue curve). Moreover, this difference between both conical 
shadow models does not increase linearly with time. By way of conclusion, 
debris trajectories obtained with our symplectic integrator coupled to our 
smooth shadow function turn out to be consistent with the ones computed 
by NIMASTEP (ABMIO) using um- 

Eventually, a last remark is given about the energy conservation in the case 
of cylindrical and conical shadow models. The comments about the quasi- 
symplecticity are, of course, the same for the solar radiation pressure than for 
the third body contributions. Hence, we are still limited in the computation 
of the relative variation of the energy. However, it is shown in Fig 11 that 
the relative error in energy does not increase with time, even on an extremely 
long time span (5000 years). Let us remark that the energy computed takes 
into account the contribution of the solar radiation pressure with permanent 
sunlight. Indeed, the part of the Hamiltonian function corresponding to the 
solar radiation pressure cannot be retrieved from the equations of motion 
i^u,p(i') Vr'Hsrp(r, 9) , the latter being impossible to integrate analytically. A 
close look to the relative error in energy shows small perturbations during 
each shadow season but it does not result in a long term drift on the energy. 
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Figure 9: Absolute difference between the orbital elements of space debris subject to the 
Earth's central attraction and solar radiation pressure with i/c and i^p shadow functions. 
The initial semi-major axes are 42164 km, the other elements are set equal to zero and 
the area-to-mass ratio is equal to 20 m^/kg. This figure emphasizes the beginning of the 
difference between both orbits after the first shadow season represented by the gray zone. 
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Figure 10: Absolute difference between the orbital elements of space debris subject to the 
Earth's central attraction and solar radiation pressure with respectively Vp and Vc shadow 
functions (in blue) and fp and j/m functions (in red) . The orbit with f m has been computed 
by NIMASTEP (ABMIO). The initial semi-major axes are 42164 km, the other elements 
are set equal to zero and the area-to-mass ratio is equal to 20 m^/kg. Each numerical 
integration has been performed with time steps equal to 150 s. 




Time [yr] 

Figure 11: Evolution of the relative error in energy of space debris motion subject to solar 
radiation pressure. The initial semi-major axis is set at 42164 km, the other elements are 
set equal to zero and the area-to-mass ratio is equal to 0.1 m^/kg. The shadow function 
fp is used to model the Earth's shadowing effects. 

5. Conclusion 

An efficient symplectic integration scheme has been presented in order to 
compute space debris motion. The underlying algorithm has been described 
and the accuracy of the integrator has been demonstrated by means of nu- 
merical comparisons. It has been pointed out that large time steps could 
be used and that the relative error in energy was not increasing with time, 
even on huge time spans. It should also be mentioned that our method is 
not stuck to one particular order and that it can adapted to the complexity 
of the perturbations and to the desired precision. Our integration scheme 
is able to take into account the Earth's gravitational potential, luni-solar 
and planetary gravitational perturbations and solar radiation pressure. In 
conclusion, it turns out that our algorithm represents a fast and reliable al- 
ternative to compute space debris or artificial satellites orbits, especially on 
long time scales. 

We have also described an innovative method to model both cylindrical 
and conical Earth's shadow crossings by means of smooth shadow functions. 
We have explained why these ones were particulary convenient in the frame- 
work of symplectic integration. It has been shown that the cylindrical model 
was not a suitable approximation of conical shadows, especially in the case 
of space debris associated to high area-to-mass ratios. It has been noticed 
that one drawback of the computation of the shadow functions during the 
integration process relies in the obligation to reduce step sizes. Even if it can 
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be considered as a limitation, one has to keep in mind that the order of the 
integrator can be reduced in order to speed up the calculations. 

Future work will be devoted to stability studies involving shadowing ef- 
fects. Such a task will be made easier by the fact that our shadow function 
derivatives are also smooth, enabling us to compute direclty the solutions of 
the variational equations. These relations are necessary to use chaos indica- 
tors like the MEGNO (Cincotta et al., 2003). Hence we will be able to extend 
the work realized in Valk et al. (2009) by considering Earth's shadows. 
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Appendix A. Orbital entrance and exit from the Earth's shadow 

In Escobal (1976), the entrance and exit true anomalies from the Earth's 
shadows are found to be the non-spurious solutions of the following function 

i?^(l + ecos/)^+/(/3cos/ + ^sin/)2 = (A.l) 

where / is the true anomaly, p = a(l — e^) is the semi-latus rectum and /3 and 
^ depend on the geocentric Cartesian position of the Sun and the inclination 
and longitude of the ascending node of the debris. As explained in Escobal 
(1976), relation (A.l) corresponds to a quartic polynom in the cosine of the 
true anomaly. Hence, this function is transformed to standard form to find a 
new quartic polynom in / which is solved in closed form by quadratic radicals 
(see Descartes' rule presented in Escobal (1976)). 

The introduction of equinoctial elements ke = ecos(f2 -|- u) and he = 
esin(f2 + u) has been proposed in Valk & Lemaitre (2008) to avoid singular 
orbital elements. It yields the following function 

Rl{l + ke COS f + he sin ff + p^{(3 cos / + ^ sin fY - = (A.2) 

whose solutions are found using the so-called resultant method (see Gronchi 
2005 for an application example) . The latter lets us solve analytically the 
problem as a system of two algebraic equations in two variables. 
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We hereby propose a further improvement in the computation of the so- 
lutions of (A. 2). By defining T := tan(//2), (A. 2) can be written as 



Rl 1 + h 



1 + T2 



T2 , 2T 

— + hez 



2^2 



rp2 2T 



2^2 



J^2 



p 



or, equivalently, 



2d2] 



+ T^[-AP^p^ + 4KRl-4KKRl] 
+ T2[4eV - 2p2/32 _ 2p2 + 4/^2/?^ + 2Rl - 2klR, 
+ T[4/3ep' + 4/ie/?^ + 4/iAi?|,] 
+ - p2 + ^1 + 2A;ei?| + klRl] 

= 

Descartes' rule can then be used to find orbital entrance and exit true 
anomalies. The advantage of this method is twofold. First, the use of the 
tangent function directly indicates the right quadrant for the true anomaly. 
Then, the resultant method is not necessary anymore. 

Appendix B. Computation of A.<p 

Relation (8) tells us that 



cos (f)2 + COS a 



Rn 



cos^ a + 



Rn 



sm a 



(B.l) 



cos (ps + COS (3 



— ^ cos^ p sm p 



0. 



(B.2) 



Then, we denote by p the small quantity r/r©. It leads to the following 
simplified expression 



r — 



01 



1 — 2p cos + ^ 1 — 2p cos I 
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where the small term has been neglected. Further calculations lead to the 
following expressions, keeping only terms of order p : 



sin a ~ tan a p ■ 



sin (3 ~ tan /3 — p 



Rq + Rq 



and cos a ~ 1 



and cos /3 ~ 1. 



Replacing (B.3) and (B.4) values in (B.l) and (B.2), one obtains 

Rq — Rq 



COS</)2 



l-^ + pR, 



cos</)3 + a/1 
By (B.5) and (B.6), we get 



—r, — P-^e 5 — 




0. 



(B.3) 
(B.4) 

(B.5) 
(B.6) 



J JO R®R& 
COS 03 - COS 02 = 2p — . 



Eventually, it can be shown that 



()2 - 03 = + O(p(0 



(B.7) 



where — 02 is small at the shadow entrance. 

Appendix C. Computation of Ah. 

At the cylindrical shadow entrance, the angle is equal to 0i. From (6) 
and (7), one obtains 



■5u(0l) = T COS 01 + COS a 



Sp(0i) = r COS 01 + COS P 



Ri, cos^ a + Rq si 



sma 



— Rq cos^ /3 — Rq sin /3 



Then, from (B.3) and (B.4), it follows that 



Suih) = r cos 01 + \/r'^ - Rq + RqP 
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Rq — Rq 



0{p' 



sp(0i) = r cos 01 + ^r^ - Rl - R^p ?^l±3^ + o{p^). 
Eventually, (8) yields 

s.{<P^) = R^p^^^^^ + Oip") 
r 

sM = -R^p^^^^±^ + 0{p'). 
r 

In conclusion, we can write 

R 

Ah = Su(0i) - Sp(0i) ~ 2pR(s — . 

r 

Appendix D. Computation of At in the non-circular case 

The only difference with the circular case is given by the dependence of r 
on (j). Expressing r in terms of Keplerian elements, one has 

ail - e') 

r{(f)) 



1 + ecos /(</)) 

where / is the true anomaly. Expanding this relation (see e.g. Murray & 
Dermott 1999) and keeping only terms of first order in eccentricity, r can be 
finally written as 

r(0) ~ a(l — e cos M(0)) = a(l — e cos(0 + ^/jq)) 

where M is the mean anomaly and tpo is the appropriate phasing. 
Following the same scheme as for the circular case, we end up with 

At = ^^(1 + e cos(0 + ^o)) + 0{p^e + pe^ + p^e^) days 



and 

Ah = -^-^-"^ [i + e cos(0 + tPo)) + 0{p'e + pe' + p'e 



2pRqRq /j, I ; \\ I /n^ 2 I 2 I 2 2\ 
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